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tl.  abstract 

This  grant  supported  research  by  two  University  of  California  Berkeley 
graduate  students  on  Innovative  numerical  methods  for  flows  around  complex 
boundaries  such  as  occur  In  solldiflcaclon  from  the  melt.  New  two-dimensional 
magnetlzatlon-based  methods  for  complex  fluid  flows  were  developed  and 
analyzed.  New  adaptive  finite  difference  methods  for  simulation  of  the 
fluid  flow  In  Czochralskl  growth  were  implemented  and  tested. 
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Thb  jjraiil  »upporte<|  research  by  two  UC  Berkeley  graduate  students  on 
computational  fluid  dynamics  and  materials  science.  Bolli  arc  developing  Innc^ 
vative  numerical  methods  for  flows  around  compWx  boundariee  such  ao  occur  in 
solidification  fiom  the  uicli.  Research  summaries  for  each  student  follow. 

RicAtdo  Cortes  (currently  an  NSF  postdoctoral  fellow  at  Courant  Insti¬ 
tute)  studied  magnetization  methods  for  fluid  flow. 

A  new  formulation  of  S-diiiiensional  Incompressible  fluid  flow  has  been  pre¬ 
sented  recently  with  the  introduction  of  a  new  variable  aometimes  lefeired  to 
«  "velocity  magneiUation”  or  simply  •magnetization  A  dUcretiiation  of  the 
resulting  equations  produces  a  Ilaiiiilioiilaii  system  which,  In  turn,  leads  to  a 
new  family  of  r.agrangian  numerical  methods  for  the  study  of  incompressible 
flow.  These  numerical  methods  ate  essentially  independent  of  upatial  dimension 
and  have  the  property  thot  they  preserve  disciete  invarianls  associated  with 
the  Hamiltonifm,  llowever  them  xre  many  issues  related  to  this  new  formula- 
lion  that  require  further  investigation  in  order  to  arrive  at  efficient  numerical 
algorithms  that  may  be  used  in  &  variety  of  situations.  These  Issues  include: 

1*  The  design  oj  robust  hybrid  tnfihni$  yiiug  eorficify  end  magnctirafion 
variahhs.  The  lUHguetization  variables  can  be  Used  in  cordunction  with 
vortex  methods  in  order  to  exploit  the  advantages  that  cocli  representation 
offers.  Accurate  and  relatively  inexpensive  procedures  to  switch  from  one 
computational  variable  to  the  other  are  re<|Ulred. 

2.  Froetdurts  in  ensiirt  that  iht  magnetigation  remains  optima/  or  nearly 
optimal.  For  a  given  initial  configuration,  magnetization  naturally  tends  to 
evolve  into  allowable  configurations  that  depart  from  the  Optimal  one.  In 
three  dimensions,  finding  the  latter  is  a  prohlem  related  to  finding  minimal 
surfaces  attached  to  a  fixed  frame.  The  accuracy  of  the  numerical  method 
is  related  to  the  distribution  of  magnetizatlou,  making  the  problem  of 
finding  the  optimal  configuration  a  crucial  pari  of  the  method. 
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8.  A^^j^licahUiiy  ©/  ihtic  mcihods  near  hoitndarics.  The  formulation  of  the 
NavIftr-StnVfts  equations  in  terms  of  inaj^uatlsation  hac  the  potential  of 
suggesting  new  ways  of  dealing  with  creation  of  vorticity  inside  boundary 
layers. 

4.  D^^inyularixalion  of  the  equations  of  motion  end  aifproyrUic  foot  a«mma« 
lion  Xtchntqnts,  The  differential  equations  that  govern  the  evolution  of 
the  particle  positions  and  the  magnetization  they  carry  have  high-order 
jvingn  lari  ties,  0(r"+*)  whore  tt  is  the  dimension  of  the  problem.  The  or¬ 
der  of  the  singularities  has  an  adverse  effect  on  the  speed  up  attained  by 
some  foet  euminotion  techniques;  therefore,  ofi  effoit  must  be  made  to 
reduce  the  order  of  the  singularities  in  the  equdtioqs  and  to  develop  fast 
sximmatioii  pioccdiircs  that  arc  well  suited  for  these  equations. 

5.  Maihemeiical  theory.  The  relationship  between  parameters  of  the  itielliod 
and  other  necessary  conditions  for  the  convargaurp  of  the  discrete  variables 
to  their  continuous  counterparts  are  only  conjectured.  There  is  a  need  for 
theory  that  provides  convergence  and  ctalility  proofs  and  error  bounds 
based  on  discretization  parameters. 

During  the  academic  years  ,191)3-114  and  1994-95,  this  AASFRT  grant  sup¬ 
ported  intensive  work  on  some  of  the  above  issues,  coiiUmiing  the  work  on  vertex 
methods  whirh  wa«  done  under  the  support  of  AFOSR  Grant  No.  FDF40620- 
93-1-0058  during  the  last  few  years.  A  summary  of  Dr.  Cortez's  results  follows. 

In  two  dimencions,  the  relationship  between  discrete  maguetizatioa  variables 
and  vortices  has  been  established  by  %  rigorous  interpretation  of  magnetization 
as  vorCcA  dipoles  with  a  prescribed  dipole  moment.  For  a  given  discretization, 
each  magnetization  vector  can  be  replaced  by  a  pair  of  vortices  of  equal  but 
opposite  strength  located  some  small  di.^ttancA  apart.  Then,  vortices  sufficiently 
•near**  others  can  be  cuiubiiicd  into  a  single  vortex  with  strength  equal  to  the 
sum  of  individual  strengths.  This  collapse  of  vortices  keeps  the  computation 
feasible.  Error  analysis  provides  a  criterion  for  how  near  vortices  should  be  in 
order  to  be  combined.  The  distance  between  the  vortices  is  a  parameter  which 
i«  chosen  so  that  the  error  in  this  procedure  io  of  a  prc4?cribed  size. 


Example  of  eombinatiou  of  vortices 

The  reverse  procedure  hoc  alco  been  accomplished.  Given  a  set  of  vortices 
of  net  vorticity  zero,  one  can  extract  subseta  to  be  represented  by  a  string 
of  vortex  pairs,  whicli  ate  then  aiiproxiinated  by  magnetization  vectors.  The 
errors  in  both  procedures  arc  shown  to  be  0(/4*),  where  h  is  the  dipole  distance. 
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The  ability  to  cwitch  from  mngncfti^at.ion  to  vorlicity  and  back  allows  iiimietical 
experimcntfi  that  compare,  particle  positions  obtained  by  magnetization  methodc 
with  those  obtained  by  well  understood  vortex  methods^  and  provides  a  very 
important  first  step  toward  the  making  of  hybrid  methods.  Tlie  inuacrical 
experimenlB  show  the  order  of  convergence. 

In  many  pliyaical  altuatioiis,  a  magnetization  couflguration,  whicll  initially 
reptosonU  a  given  set  of  vortices  to  a  preccribed  accuracy,  evolves  naturally  in  a 
way  that  endures  accuracy  loss  to  the  point  where  the  magnetization  no  longer 
opproximates  the  original  voiticcs  evolved  111  time.  The  velocity  field  due  to 
a  string  of  magnetization  vectors  connecting  two  vorticco  has  been  interpreted 
as  an  approximation  to  the  line  integral  that  represents  the  velocity  field  <!iie 
to  the  vorticcfi.  The  low  of  accuracy  is  now  understood  as  the  deterioration 
of  the  approximate  line  Integral  due  to  the  stretching  of  the  curve  along  which 
wc  integrate.  This  understanding  has  led  To  various  procedures  that  correct 
the  problem  and  encure  that  the  magnetization  remains  optiiiiAl  by  inaiiitaining 
the  initial  accuracy  over  time.  Some  of  thase  procedures  periodically  discard 
the  integration  curve  and  replace  it  with  a  more  convenient  one,  computing  a 
new  magnetization  field  along  the  new  curve.  Others  use  the  original  cuivc  but 
refine  the  intervals  used  in  the  quadrature  so  aa  to  maintain  the  initial  accuracy. 
The  Utter  case  is  leas  efncieiit  , since  the  number  of  particles  may  grow  without 
hound,  hut  it  has  a  natural  extension  to  3-dimcncioiial  problems.  There,  a 
magnetization  vector  is  approximated  by  a  vortex  loop  on  the  plane  norma]  to 
the  vector^  thcac  loops  can  be  split  into  smaller  uiies  to  iiialul.ain  the  accuracy 
of  the  discretization.  The  first  draft  of  the  work  above  has  been  written  os  port 
of  Mr.  Cortez's  dissertation. 

Magnotizotion  variables  are  particularly  well  auiteil  for  problems  in  which 
vort-cx  dipoles  are  used.  One  such  problem  is  the  motion  of  an  clastic  membrane 
iiaiiiei^  ill  a  fluid,  where  the  effects  of  elastic  forces  acting  on  the  fluid  can 
he  introduced  the  rate  of  change  of  magnetization.  The  forces  that  the 
membrane  imparts  on  the  fluid  over  a  small  time  inUrval  represent  impulse, 
and  their  cffecl  can  be  introduced  via  vortex  dipoles.  In  other  words,  the  forces 
represent  an  evolving  magnetization  field  on  the  membrane.  This  applicolloa 
has  been  implemented  in  two  dimensions  for  smooth  membranes  and  invUcid, 
incompressible  flow. 

is  thought  of  os  a  starting  point  for  the  understanding  of  magnetiza¬ 
tion  and  for  providing  extensions  to  K^.  In  three  dimensions,  magnetization 
variables  have  the  advantage  over  vortex  methods  of  automatically  yielding 
a  divergence-free  vortlcity  field.  This  ie  a  source  of  problems  with  3-D  vor¬ 
tex  methods.  Problems  with  3-D  magnetization  metliod.s  include  the  arbitrary 
shape  of  the  loops  when  two  neighboring  loops  must  be  combined,  and  the  need 
to  find  minimal  surfaces  when  a  new  magnetization  rcprccciitotlon  of  large  vor¬ 
tex  loops  is  needed.  Computationally  efficient  extensions  of  the  ^dimensional 
results  will  provide  a  promising  arena  foi  tlie  luodellug  of  3-dimensional  fluid 
flow. 
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Haas  JoliAAson  U  a  doctoral  candidate  in  tltc  Department  of  Mechanical 
En^^ineerihg  at  the  University  of  California,  Berkeley.  His  r<isearc}>  fornses  on 
finite  difference  methods  for  .sininlation  of  unsteady  fluid  dyiiaiaics,  especially  In 
hulk  C^oehraUki  crystal  growth  of  multi-compon^nt  semiconductor  eubatralcs. 
One  of  the  major  factors  that  limits  the  quality  of  such  crystals,  is  the  effect  of 
convection  in  the  melt  on  thermal  gradients  and  composition  in  the  solicl  phaae 
(for  example,  see  Derby  (2J).  Ry  armrately  representing  the  melt  int^rfaca, 
fluid  dynamics,  and  heat  transfer,  one  can  belter  predict  ctysr.al  cliaracteristics 
based  on  the  furnace  geometry  and  thermal  environment.  In  addition,  because 
the  growth  process  is  lengthy,  there  is  the  possibility  of  artive  control  of  certain 
crucible  factors,  sucb  u  ai>])Iicd  iiiagiietic  fields  and  time-dcpeiident  heating. 

Mr.  Johansen*a  approach  U  based  on  previous  work  with  ‘‘Cartesian  grid” 
and  *volume-of-fliud"  methods  for  unsteady  fluid  dynamic.  These  methods 
represent  the  pioMem  on  a  square  finite  difference  grid  and  apply  standard  dis¬ 
cretization  techniques  to  the  interior  of  the  domain.  Near  irregular  boundaries, 
on  a  smaller  number  of  points,  special  algorithms  are  used  to  advance  the  solu¬ 
tion,  taking  the  local  geometry  into  account,  Tlitajc  Cartesian  grid  luclliods  were 
recently  applied  to  invisrid  compressible  flows  by  Pembet,  «t  al.[3].  Time  step 
constraints  (due  to  small  cells)  were  overcome  by  redistributing  excess  fluxes  to 
larger  neighboring  cells,  while  also  maintaining  conservation.  The  appiuacli 
also  includes  AMR,  local  Adaptive  Mesh  Refinement  of  the  finite-difference 
grid,  bal^ed  Oil  work  done  by  Berger  and  Colella  [4];  results  for  two-  and  tliree- 
dimeneiona]  regions  compare  well  with  other,  more  expensive  computations.  A 
similar  approach  has  been  developed  for  {nromptessihle  flows  by  Almgren  «t 
al.  (13],  wlilcli  included  efficient  solution  of  a  Poisson  equation,  with  no-flux 
boundary  conditions,  on  the  irregular  domain,  Volumc-of-fluid  algorithms  have 
been  applied  ic  a  number  of  problems  in  compressible  flow.  Chern  and  Colella 
[6]  implemented  conservative  iuteifocc  tracking  for  compressible  flow,  using  a 
Simple  Line  Interface  Calculation  (SLIC),  to  represent  the  boundary  between 
different  gases.  When  combined  with  AMR,  the  method  produced  excellent  cor¬ 
relation  with  experiments  of  shocks  impinging  on  gas  tut'er faces  [C],  even  with 
the  first-order  representation  of  the  front.  Since  then,  cimilat  volume  of  fluid 
algorithms  have  been  developed  for  multifluid  compressible  flows  [7],  even  in¬ 
cluding  combustion  [8],  and  complicated  physical  domains  [9]. 

Thei^e  papers  demonstrate  that  modeling  the  melt  motion  is  certainly  within 
the  capabilities  of  a  volume-of-fluid  method.  The  largest  hurdle  in  applying  it 
to  cryetal  growth  is  developing  an  algorithm  for  tracking  the  iiiteiface,  while 
maintaining  con.servation  of  energy.  Once  this  is  accomplished,  a  multl-etep, 
prcilictor-corrector  algorithm  can  be  used  to  advance  the  solution  in  time,  in 
which  both  the  melt  interface  and  container  geometry  will  be  represented  using  a 
hybrid  front-tracking  method  derived  from  work  mentioned  above.  Combining 
the  Stefan  problem  with  an  integrator  for  the  Boussiaesq  equations  for  the 
melt  (  similar  to  recent  progress  for  Incompressible  flows  [11,  13,  14]  ),  will  be 
the  last  step  for  the  base  algorithm.  Adding  an  adaptive,  finite-difference  grid 
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hietar<;liy  will  provide  locM  refinement  of  iateifdcett^  and  lidp  verify  the  method. 
Implementation  of  tbo  rornplex  dgorithm  Is  simplified  with  the  uco  of  hybrid 
C++  and  Forlraji  libraries  developed  at  Lawrence  Livermore  National  Labs. 

The  details  of  the  method  arc  hoacd  on  tlie  papers  nicntioncil  above,  with 
adaptations  and  improvements  for  the  St4»fan  problem.  Our  front  recon8true> 
tion  algoiithm  will  be  hiujcd  on  work  done  by  Pilliod  and  Puckett  [10].  Their 
approach  provides  a  method  of  approximating  an  interface,  based  on  a  linear, 
least-squares  fit  to  the  local  array  of  volume  fractions.  The  algorithm  is  shown 
to  be  second-order  accuiate  In  »r»acc,  even  when  applied  to  complicated  shapes. 
For  the  Stefan  problem,  the  iuterfice  mov«c  normal  to  itself,  with  a  speed  pro¬ 
portional  to  the  local  temperature  gradient.  Chern  and  (Tolella  [fi]  implemented 
ohoclc  tracking  for  compressible  flow,  which  is  directly  applicable  to  the  motion 
of  the  front  in  the  Stefan  problem.  By  coupling  the  reconetruction  algorithm 
with  an  operatot-splli  version  of  ibis  front-tracking  scheme  (  modified  to  he 
minimum  and  maximum  preserving  ),  we  will  be  able  to  advance  the  phase 
boiindary  in  time,  f  Conservation  of  energy  is  enforced  by  adding  equivalent  heat 
sources,  due  to  ciitlialpy  of  formation,  Where  necessary. 

Solving  the  temperature  equation  on  the  irregular,  changing  domain  will  be 
the  next  step.  Ir'irst,  the  lialf-time  velocitie.5t  will  be  found  using  an  algorithm 
similar  to  those  in  [13].  An  in(«rmediate  MAC  projection  guarantees  that  the 
edge-centered  velocities  have  aero  diectete  divergence,  making  them  ideal  for 
computing  a  time-centered »  temperature  convection  term,  l^iis  will  then  be 
used  in  the  right-haud  side  of  the  temperature  equation,  along  with  the  en¬ 
thalpy  of  formation  mentioned  earlier.  The  melting  point  of  the  substance  wHl 
require  a  Dirichict  boundary  condition  for  temperature  on  the  interface;  how¬ 
ever,  prescribing  the  front  motion  is  ctpiivalcut  to  specifying  the  temperature 
gradient  at  the  interface,  too.  This  conflict  will  be  surmounted  by  ignoring  the 
constraint  on  the  temperature  gradient;  it  Is  determined  by  fitting  a  quadratic 
through  values  interpolated  from  full  cells  near  the  lii1.eifacc,  using  the  normal 
and  midpoint  of  the  reconstructed  front.  The  implicit  discretlration  of  the  difTii- 
sion  equation  will  be  solved  using  mtiUi-grid  iterations  to  accelerate  convergence 
[11]  [13]  .  Once  the  new  temperature  field  is  determined,  the  energy  loss  due 
to  a  niismatcli  in  gradients  will  he  added  into  the  temperature  equation  on  the 
next  tune  step. 

Finally,  the  front  position  and  temperature  fields  will  be  used  in  the  inlo- 
gration  of  the  Boussinesq  equationa  for  flow  in  the  melt.  The  remainder  of  the 
algorithm  is  similar  to  tliat  presented  in  (13l  and  (14).  The  main  differences 
will  be  the  inelusion  of  vUcosily,  which  demands  no-slip  boundary  conditions  at 
all  interfaces.  Tlie  viscous  forces  will  be  treated  similarly  to  the  heat  equation, 
with  the  advection  terms  and  body  forces  (due  to  gravity  and  externally  applied 
magnetic  fields)  serving  os  sources  in  the  Crank  Nicholson -type  discretization. 
In  the  final  step,  the  velocity  field  will  be  projected  onto  its  divergence-free  part. 

Progrcaa  to  date  has  foc«.s6cd  on  acciirate  integration  of  the  Stefan  prol>- 
lem,  and  the  heat  equation  with  Dirichlst-type  boundary  conditions.  First, 
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Mr*  Joh^l'iSen  has  thoroughly  Usted  i.  one>dini«n&ion&I  method  for  the  cl&sslcat 
Stefan  iHoldrin,  which  reprcseriu  a  balance  between  heat  conduction  and  en¬ 
thalpy  of  formation.  The  interface  gradient  is  calculated  using  a  quadratic  fit 
to  nearby  temperature  values;  after  the  froiit  is  advanced,  conservation  is  en¬ 
forced  by  adding  apiiropriate  source  terms*  The  quadratic  temperature  fit  and 
lagged  conservation  have  been  shown  to  be  stable  for  ctiiaiu  diMTetizatlons; 
results  dernonstrat<e  second-order  convergence  in  space,  and  first-order  accuracy 
in  time,  regardless  of  Stefan  mimber.  Second,  the  solution  of  Diridilet-lype 
Poisson  problems  on  irregular,  two>^d{nioncioiial  domalne,  provides  part  of  the 
tirne^stepping  algorithm  for  the  Stefan  problem  in  more  than  one  dimoueion. 
The  calculation  of  interface  gradients,  luujcd  on  quadratic  fits  to  interior  val¬ 
ues,  has  been  demonstrated  as  globally  second  order  accurate,  even  for  front 
lepicsrutatious  containing  arbitrarily  small  cells, 

Mr.  Jobaneau  will  have  a  working  code  by  December,  1990,  Adaptive  so¬ 
lution  of  the  Stefan  problem  in  multiple  dimensions  will  be  fininlicd  by  June, 
1996;  it  will  piovldc  the  crucial  step  needed  to  incorporate  the  front  motion  into 
a  Cartesian  grid,  Incoinprecciblc  Navicr-Stokes  solver,  whlih  is  being  developed 
at  both  UC  Berkeley  and  LLNL,  and  will  be  completed  by  September,  1996. 
The  remaining  three  months  will  be  spent  lasting  the  algorithm,  by  verifying 
it  against  experimental  and  computational  rcsulta;  cxtrrnally-applie<l  magnetic 
fields  might  be  included  in  the  model  in  Spring  of  1997.  This  work  will  be  a 
major  step  in  applying  finitc-difieieiice  algorithms  to  real  engineering  problems. 
The  simulation  of  semiconductor  substrate  crystal  growth  will  be  a  new  appli- 
cation  of  the  Cartesian  grid  method,  with  front  tracking,  conv<>rtion  of  heat  and 
melt  compoGition,  and  adaptive  icfinetiieiii  of  the  finite-difTcrcnce  grid  around 
significant  features. 
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